
%***********************************************************************
% run_struc_model.m
%***********************************************************************
% C.-Y. Cynthia Lin, 27 July 2011
%-----------------------------------------------------------------------
% The code provided is from the following paper:
%
% Lin, C.-Y. Cynthia.  (2013).  Strategic decision-making with information 
% and extraction externalities:  A structural model of the multi-stage 
% investment timing game in offshore petroleum production.  Review of 
% Economics and Statistics, 95 (5): 1601-1621.
%
% If the code is helpful to you in writing your own code, I would appreciate
% it if you would cite the above paper in your work and acknowledge me 
% for help with the code.  Thank you!
% ------------------------------------------------------------------------
%
% This code estimates the structural econometric model in the following 
% paper:
%
% Lin, C.-Y. Cynthia.  (2013).  Strategic decision-making with information 
% and extraction externalities:  A structural model of the multi-stage 
% investment timing game in offshore petroleum production.  Review of 
% Economics and Statistics, 95 (5): 1601-1621.
%
% 9 parameters are estimated (i.e., profit fxn is linear and there are
% no interaction terms), including coefficient on drilling cost
% in exploration profit function
%
%
% STATE VARIABLES
% ---------------
% The state variables are:
%	tot_e_t = total number who explored before lease-yr t
%	tot_d_t = total number who developed before lease-yr t
%	bid_bin_t = discretized avg high bid (avged over the 2 tracts)
%	ce_bin_t = discretized real drilling cost 
%	oilp_bin_t = discretized real oil price
%
% DATA
% ----
% The variables in the panel data input file (created by the
% idl program make_panel.pro) are:
%	mkt_i 
%	t		= yr since lease sale
% 	state variables:
%		tot_e_t = total number who explored before lease-yr t
%		tot_d_t = total number who developed before lease-yr t
%		bid_bin_t = discretized avg high bid (avged over the 2 tracts)
%		ce_bin_t = discretized real drilling cost 
%		oilp_bin_t = discretized real oil price
%       tot_not_e_t = total number who have not explored before lease-yr t 
%       tot_e_not_d_t = total number who are explored at t (ie, by t+1)
%           but have not developed before lease-yr t 
%       n_e_t 		= number who explore from t to (t+1)
%	n_not_e_t 	= number of unexplored who do not explore 
%				from t to (t+1) 
%       n_d_t 		= number who develop from t to (t+1)
%       n_e_not_d_t 	= number who are explored at t (ie, by t+1)
%				but undeveloped before t who 
%				do not develop from t to (t+1)
%       n_e_and_d_t 	= number who both explore & develop from t to (t+1)
%
% The variables we will use from the profit data input file 
% (created by the idl program make_panel.pro) are:
%	gross_profits = (1/inv_royalty_rate) * revenue - cost
%			   for those tracts that develop
%
%
% PARAMETERS
% ----------
% The parameters to be estimated are:
%     sigma_mu
%     sigma_eps
%     ce_factor
%     gamma_tote
%     gamma_totd
%     gamma_bid_bin
%     gamma_cd_bin  (assume ce_bin = cd_bin)
%     gamma_oilp_bin
%     gamma_constant
%
%
%
% MOMENT CONDITIONS
% -----------------
% The moment function for each market-time observation j is: 
%    f(data_j,params) = ...
%      [ ge(tuple_j, params) - ge_emp_j(tuple_j) 
%	 gd(tuple_j, params) - gd_emp_j(tuple_j) 
%	omega_j * ( ge(tuple_j, params) - ge_emp_j(tuple_j) )
%	omega_j * ( gd(tuple_j, params) - gd_emp_j(tuple_j) )  ]
%
%
% The population moment conditions are given by:
%     E[f(data_j,params)]=0 
%
% The sample moment conditions are given by:
%    (1/n sum_j{f(data_j,params)})=0 
% and amount to:
%    sample_mom_conds = 
%    (1/n sum_j{f(data_j,params)})=
%      1/n *
%       [ sum_s{tot_not_e_tuple(s) * (ge(s) - ge_emp(s))} 
%         sum_s{tot_e_not_d_tuple(s) * (gd(s) - gd_emp(s))} 
%         etc     ]
% where s indexes the tuples reached in the data
%
% The weighted moment conditions is (1 X 1)
%    (sample_mom_conds)' weight_matrix^(-1) (sample_mom_conds]) 
%
%
%***********************************************************************


% Variables to be modified by user

    % Boolean for whether want tracts in market to be owned
    % by the same firm
    use_same_firm = 0;

    % Acreage of tracts to use:
    acreage_type = 0; % all acres
    

% Constants

    % Lease term
    lease_term = 5;

    % Discount factor
    beta = 0.9;

    % Number k of parameters
    n_params = 9; 

    % # bins for bid_bin 
    n_bins_b = 3;

    % # bins for ce_bin
    n_bins_c = 2;

    % # bins for oilp_bin
    n_bins_p = 3;

    % Number of iterations in the optimization (both coarse & fine)
    maxit_coarse = 1000;
    maxit_fine = 1000;
    
    % Termination tolerance on the moments in determining 
    % convergence (both coarse & fine)
    tol_moms_coarse = 5; 
    tol_moms_fine = 1e-2;  
    
    % Termination tolerance on the parameters in determining 
    % convergence (both coarse & fine)
    tol_params_coarse = 5; 
    tol_params_fine = 1e-2; 

    % Termination tolerance when determining the fixed point
    % for the value functions (both coarse & fine)
    tol_fp_coarse = 20; 
    tol_fp_fine = 1e-2; 


    % Initial difference between value functions before solve
    % for fixed point
    ini_diff_fp_coarse = 1.5*tol_fp_coarse;
    ini_diff_fp_fine = 1.5*tol_fp_fine;

    % Vector of Vce estimators
    Vce_estimators_vec = cellstr(['G']);
    

    % Initial guesses of parameters
    sigma_mu_ini		= 5;
    sigma_eps_ini		= 5;
    ce_factor_ini		= -10;
    gamma_tote_ini		= 0;
    gamma_totd_ini		= 0;
    gamma_bid_bin_ini		= 5;
    gamma_cd_bin_ini 		= -10;
    gamma_oilp_bin_ini		= 5;
    gamma_constant_ini		= 5;
	  

    
    % Last period of consideration for the exploration stage
    % (this is the last period in which decisions can be made)
    Te = lease_term - 1;

    % Last period of consideration for the development stage
    % (all years from lease_term onwards are lumped together
    % b/c tot_e_t, which depends on t, will not change from then
    % onwards)
    Td = lease_term;


% Determine samefirm_str
switch use_same_firm
  case 1
    samefirm_str = '_samefirm';	  
  otherwise
    samefirm_str = '';	  
end

% Determine acreage-type
switch acreage_type
  case 1 
    acreage_type_str = '_5760acres';
  case 2 
    acreage_type_str = '_5000acres';
  case 3 
    acreage_type_str = '_lt5000acres';
  case 4 
    acreage_type_str = '_ge5000acres';
  case 5 
    acreage_type_str = '_lt4000acres';
  case 6 
    acreage_type_str = '_ge4000acres';
  case 7 
    acreage_type_str = '_lt3000acres';
  case 8 
    acreage_type_str = '_ge3000acres';
  case 9 
    acreage_type_str = '_4995acres';
  case 10
    acreage_type_str = '_2500acres';
  otherwise
    acreage_type_str = '';
end


% Make vector of initial guesses for the parameter values
params0 = zeros(n_params,1);
params0(1,1) 	= sigma_mu_ini;
params0(2,1) 	= sigma_eps_ini;
params0(3,1) 	= ce_factor_ini;
params0(4,1) 	= gamma_tote_ini;
params0(5,1) 	= gamma_totd_ini;
params0(6,1) 	= gamma_bid_bin_ini;
params0(7,1) 	= gamma_cd_bin_ini;
params0(8,1) 	= gamma_oilp_bin_ini;
params0(9,1) 	= gamma_constant_ini;

       
% Directory and file names
    % Directory for programs
    programs_dir = ...
      strcat('C:/CYCLin_files/2002-2003/oil_investtiming/programs/', ...
             'structural_model/struc_model_v8'); 

    % Directory & file with input data 
    input_dir = 'C:/CYCLin_files/2002-2003/oil_investtiming/panel_data/'; 
    panel_file = strcat('panel_NS6_EW6_1980_2tract', ...
    			acreage_type_str, ...
                        samefirm_str, ...
			'.dat');
    profit_file = strcat('profit_NS6_EW6_1980_2tract', ...
    			acreage_type_str, ...
                        samefirm_str, ...
			'.dat');

    % Directory & file to output estimated parameters 
    output_dir = 'C:/CYCLin_files/2002-2003/oil_investtiming/results/'; 
    %output_file = 'oilmodel_v7_NS6_EW6_1980_2tract.out';
    %output_file = 'oilmodel_v7_NS6_EW6_2tract.out';

    % Directory & file to output results for use in parametric bootstrap 
    pbs_dir = ...
      'C:/CYCLin_files/2002-2003/oil_investtiming/pbootstrap_panels/'; 



% Go to input directory
cd(input_dir)    

% Open input files
panel_dat 	= load(panel_file, '-ascii');
profit_dat	= load(profit_file, '-ascii');



% Error check
[n_obs, tot_n_vars] = size(panel_dat);
disp(' n_obs'); disp(n_obs);
disp(' tot_n_vars'); disp(tot_n_vars);
disp(' n_vars'); disp(tot_n_vars - 2);

% Go back to programs directory
cd(programs_dir)    

% Convert panel_dat to panel_arr
[panel_arr, max_mkt_size] = ...
    convert_paneldat_to_panelarr(panel_dat); 

% Count various variables by tuples of state variables
[n_per_tuple_arr, ...
 tot_not_e_tuple, ...
 tot_e_not_d_tuple, ...
 n_e_tuple, ...
 n_d_tuple, ...
 n_e_and_d_tuple, ...
 n_totnote_totenotd_ne_nd_ned_tuple, ...
 n_tuple_to_tuple, ...
 n_not_e_tuple_to_tuple, ...
 n_e_not_d_tuple_to_tuple] = ...
    count_by_tuples(panel_arr, max_mkt_size, lease_term, Te, Td, ...
     	n_bins_b, n_bins_c, n_bins_p);


% Make a matrix of the state variable tuples found in the data,
% form the empirical transition matrices M 
% and the empirical probabilities g,
% and form vectors of other variables by tuple found in the data 
[tuple_mat, n_per_tuple_vec, ...
          tot_not_e_vec, tot_e_not_d_vec, n_e_vec, n_d_vec, ...
	  n_e_and_d_vec, ...
	  n_totnote_totenotd_ne_nd_ned_arr, ...
          Mn_emp, Me_emp, ge_emp, gd_emp] = ...
      form_emp_M_and_g(...
      	  lease_term, Te, Td, ...
      	  max_mkt_size, ...
          n_per_tuple_arr, ...
          tot_not_e_tuple, ...
          tot_e_not_d_tuple, ...
          n_e_tuple, ...
          n_d_tuple, ...
          n_e_and_d_tuple, ...
	  n_totnote_totenotd_ne_nd_ned_tuple, ...
          n_tuple_to_tuple, ...
          n_not_e_tuple_to_tuple, ...
          n_e_not_d_tuple_to_tuple);

% Error check
disp(' n_obs'); disp(n_obs);
disp(' n_obs from n_per_tuple_vec'); disp(sum(n_per_tuple_vec)); 


% Sum the revenues and gross profits for tracts that developed
% by tuple of state variables
[sum_revenueifd_tuple, sum_profitsifd_tuple] = ...
    sum_profitsifd_by_tuples(profit_dat, max_mkt_size, lease_term, Td, ...
    n_bins_b, n_bins_c, n_bins_p);


% Calculate the empirical average revenue and gross profit for
% tracts that develop
[avg_revenueifd_emp, avg_profitsifd_emp] = ...
  form_emp_profitsifd(...
	  tuple_mat, ...
	  sum_revenueifd_tuple, ...
	  sum_profitsifd_tuple, ...
          n_d_tuple, ...
	  Td);

% Convert the individual empirical average revenue and gross profit
% into millions of dollars
avg_revenueifd_emp = avg_revenueifd_emp / 1e6;
avg_profitsifd_emp = avg_profitsifd_emp / 1e6;


% Loop for each Vce estimator
for Vce_i = 1:length(Vce_estimators_vec)	  

    % Determine estimators
    Vce_estimator = Vce_estimators_vec(Vce_i);  
    disp(' Vce_estimator:'); disp(Vce_estimator);

    % Estimate the parameters in the model
    [sigma_mu, sigma_eps, ...
    	  ce_factor, ...
	  gamma_tote, gamma_totd, ...
	  gamma_bid_bin, gamma_cd_bin, gamma_oilp_bin, gamma_constant, ...
	  Vcn, Vce, ...
	  ge, gd, ...
	  pi0_e, pi0_d, EVn, EVe, wtd_moments] = ... 
 	GMM(lease_term, Te, Td, max_mkt_size, ...
	  beta, tuple_mat, n_per_tuple_vec, ...
          tot_not_e_vec, tot_e_not_d_vec, ...
	  n_d_vec, ...
          Mn_emp, Me_emp, ge_emp, gd_emp, ...
          avg_profitsifd_emp, ...
	  n_params, ...
	  maxit_coarse, tol_moms_coarse, tol_params_coarse, ...
	  tol_fp_coarse, ini_diff_fp_coarse, ...
	  maxit_fine, tol_moms_fine, tol_params_fine, ...
	  tol_fp_fine, ini_diff_fp_fine, ...
	  Vce_estimator, params0);


    % Output file for results
    output_file = char(strcat('oilmodel_v8_', ...
   			      int2str(n_params), 'params', ... 
                              '_Vce', Vce_estimator, ...
			      '_NS6_EW6_1980_2tract', ...
			      acreage_type_str, ...
			      samefirm_str, ...
			      '.out'));

    % Output files for params, Vcn, Vce, tuple_mat
    params_outfile = char(strcat('params', ...
                              '_oilmodel_v8_', ...
   			      int2str(n_params), 'params', ... 
                              '_Vce', Vce_estimator, ...
			      '_NS6_EW6_1980_2tract', ...
			      acreage_type_str, ...
			      samefirm_str, ...
			      '.dat'));
    Vcn_outfile = char(strcat('Vcn', ...
                              '_oilmodel_v8_', ...
   			      int2str(n_params), 'params', ... 
                              '_Vce', Vce_estimator, ...
			      '_NS6_EW6_1980_2tract', ...
			      acreage_type_str, ...
			      samefirm_str, ...
			      '.dat'));
    Vce_outfile = char(strcat('Vce', Vce_estimator, ...
                              '_oilmodel_v8_', ...
   			      int2str(n_params), 'params', ... 
                              '_Vce', Vce_estimator, ...
			      '_NS6_EW6_1980_2tract', ...
			      acreage_type_str, ...
			      samefirm_str, ...
			      '.dat'));
    tuple_mat_outfile = char(strcat('tuple_mat', ...
                              '_oilmodel_v8_', ...
   			      int2str(n_params), 'params', ... 
                              '_Vce', Vce_estimator, ...
			      '_NS6_EW6_1980_2tract', ...
			      acreage_type_str, ...
			      samefirm_str, ...
			      '.dat'));


    % Delete previous version of output (log) file;
    delete(strcat(output_dir, output_file));

    % Open output (log) file;
    diary(strcat(output_dir, output_file));

    % Display output
    Vce_estimator
    maxit_coarse
    maxit_fine
    tol_moms_coarse
    tol_moms_fine
    tol_params_coarse
    tol_params_fine
    tol_fp_coarse
    tol_fp_fine
    n_obs
    sigma_mu
    sigma_eps
    ce_factor
    gamma_tote
    gamma_totd
    gamma_bid_bin
    gamma_cd_bin
    gamma_oilp_bin
    gamma_constant


    % Close output (log) file;
    diary('off');

    % Form parameter vector
    params = zeros(n_params,1);
    params(1,1) 	= sigma_mu;
    params(2,1) 	= sigma_eps;
    params(3,1) 	= ce_factor;
    params(4,1) 	= gamma_tote;
    params(5,1) 	= gamma_totd;
    params(6,1) 	= gamma_bid_bin;
    params(7,1) 	= gamma_cd_bin;
    params(8,1) 	= gamma_oilp_bin;
    params(9,1) 	= gamma_constant;

    % Go to parametric bootstrap directory
    cd(pbs_dir)    

    % Save params, Vcn, Vce, panel_dat to respective files 
    save(params_outfile, 'params', '-ascii');
    save(Vcn_outfile, 'Vcn', '-ascii');
    save(Vce_outfile, 'Vce', '-ascii');
    save(tuple_mat_outfile, 'tuple_mat', '-ascii');
	  
    % Go back to programs directory
    cd(programs_dir)    

end
	  

